Entanglement at a Two-Dimensional Quantum Critical 
Point: a T=0 Projector Quantum Monte Carlo Study 



Stephen Inglis^ and Roger G. Melko^'^ 

^Department of Physics and Astronomy, University of Waterloo, Ontario, N2L 3G1, 

Canada 

^Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada 

Abstract. Although the leading-order scaling of entanglement entropy is non- universal 
at a quantum critical point (QCP), sub-leading scaling can contain universal behaviour. 
Such universal quantities are commonly studied in non-interacting field theories, however 
it typically requires numerical calculation to access them in interacting theories. In this 
paper, we use large-scale T = quantum Monte Carlo simulations to examine in detail 
the second Rcnyi entropy of entangled regions at the QCP in the transverse-field Ising 
model in 2-|-l space-time dimensions - a fixed point for which there is no exact result 
for the scaling of entanglement entropy. We calculate a universal coefficient of a vertex- 
induced logarithmic scaling for a polygonal entangled subregion, and compare the result 
to interacting and non-interacting theories. We also examine the shape-dependence of the 
Renyi entropy for finite-size toroidal lattices divided into two entangled cylinders by smooth 
boundaries. Remarkably, we find that the dependence on cylinder length follows a shape- 
dependent function calculated previously by Stephan et al. [New J. Phys., 15, 015004, (2013)] 
at the QCP corresponding to the 2 4- 1 dimensional quantum Lifshitz free scalar field theory. 
The quality of the fit of our data to this scaling function, as well as the apparent cutoff- 
independent coefficient that results, presents tantalizing evidence that this function may 
reflect universal behaviour across these and other very disparate QCPs in 2 -|- 1 dimensional 
systems. 
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1. Introduction 



At a quantum critical point (QCP), the Renyi [1] entanglement entropies contain the much- 
celebrated ability to access the central charge of the associated conformal field theory (CFT) 
in 1 + 1 space-time dimensions [2, 3]. This has provided the mainstay in a healthy dialog 
between field theory and numerical lattice simulations, where unbiased methods such as 
density matrix renormalization group (DMRG) [4] are able to calculate the Renyi entropies 
of order a precisely in quantum models in one spatial dimension (1-D). The form for systems 
with open (ji = 2) or periodic [rj = 1) boundaries in ID [5, 6], 



can be straightforwardly compared to lattice numerics [7], where measuring both L and x 
in terms of the lattice spacing a gives this term access to the universal central charge, c. 
The simplicity of the entangled boundary in Z) = 1 means that this and other geometrical 
factors, such as occur when Renyi indices a > 1 [8], can be compared between field theory 
and numerics to a high degree of accuracy in 1 + 1 space-time dimensions. 

A similar success is only beginning to be enjoyed in D > 1. Scalable finite-size lattice 
simulation methods which can address behaviour at critical points (where the length scale 
diverges), such as quantum Monte Carlo (QMC), predominantly measure the strongest 
signal from a looming non-universal leading-order contribution proportional to the boundary 
length [9, 10], L^-^/a^'^ = i (except in the presence of a Fermi surface, where even higher 
entanglement is expected [11, 12]). The subleading terms that exist must be obtained using 
subtraction of very large statistically-fiuctuating contributions from this "area law" (unless 
they can be calculated separately, as in infinite-lattice linked-cluster expansions [13]). What 
is left is a universal piece, which may depend on the shape or topological characteristics of the 
boundary, but is believed to not depend on the entangled volume. This geometry-dependence 
can be exploited to access different universal numbers in strongly-interacting models, but 
the types of geometries amenable to comparison between lattice-model simulations and 
continuum theories has so far been limited. 

One promising function to access universal quantities in 2+1 is the logarithmic 
contribution that comes from vertices in a polygonal-shaped entangled region [14, 15, 16, 17, 
18]. The universal coefficients, which depend on the vertex angle 6, have been compared in 
the past to calculations on finite-size lattices where the entangled region is a square, with four 
independently-contributing corners [6 = 7r/2). Particularly relevant are recent numerical 
results on interacting models [13, 19, 20], however to date, the only field theory calculations 
that have been performed are on non-interacting theories [15], inhibiting a quantitative check 
of universality. 

However, other shape-dependent contributions, that occur with smooth boundaries 
but otherwise impart important geometric characteristics of a 2+1 dimensional space-time 
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geometry, are also accessible by a simulation cell cut into two subregions (e.g. a torus 
bifurcated into two cylinders), and reveal sub- leading contributions to the Renyi entropy. 
Speculation exists about whether this shape dependence reveals an underlying universality 
[21, 22]. Quantum Monte Carlo simulations on interacting QCPs in 2 + 1 dimensions are 
uniquely poised to answer this question. 

In this paper we study the subleading shape dependence of the Renyi entanglement 
entropy on L x L toroidal lattices with spatial dimension D = 2, using the critical point of 
the transverse field Ising model (TFIM), 



where "o^j is a Pauli spin operator, accessed by a novel projector QMC algorithm that works 
strictly at T = 0. Simulating a variety of lattice sizes up to 40 x 40 reveals universal vertex 
contributions that are close to a non-interacting field theory in 2 + 1 . For smooth boundaries 
bifurcating the torus into two cylinders, subleading contributions show close functional form 
to Eq. (1) when the two cylinder lengths are associated with x and L — x. However, a cutoff- 
(or L/a-) independent coefficient is not obtained on finite lattices. Instead, if we look at 
the well-studied 2+1 dimensional quantum Lifshitz model [14, 23, 24, 25] (the field-theory 
associated with the Rokhsar-Kivelson (RK) Hamiltonian [26]) there is a proposed functional 
form [22] for the universal subleading term to the area law that is in good agreement with 
our data with a size- independent coefficient. This allows us to speculate on the universality 
of the scaling function, and the interpretation of its coefficient as a measurable witness to 
the universality class. Such potential demonstrates the importance of obtaining efficient 
simulation methods for calculation of Renyi entropies in strongly-interacting lattice models, 
and the continuing need for field theory calculations on interacting fixed points in 2+1 
dimensions. 

2. Entanglement at strongly-interacting critical points in 2+1 D 

Extensively studied in one spatial dimension, a multi-disciplinary community is beginning to 
examine the scaling of entanglement entropy in two (spatial) dimensions, thanks to the rapid 
development of both theory, and numerical methods for calculating entanglement-related 
quantities in scalable simulations. With the important exception of many-body systems 
housing a Fermi surface [11, 12], the prevailing paradigm for entanglement in ground state 
wavefunctions is the area (or boundary) law [9, 10], 



where A is a non-universal constant, i = ^/a^ ^ is cutoff-dependent, and the ellipses 
indicate a combination of different subleading corrections. Heuristically, the existence of an 
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area law can be related to the finite extent of correlations across the entangling boundary 
[27, 28]. 

In gapped systems, the subleading term may have contributions from several constants, 
some universal, some not. The most important universal subleading correction gives a 
contribution dependent on the topology of the entangled surface in a fractional topological 
phases, and is called the topological entanglement entropy: 

Sa = Ai- 7, (4) 

where e.g. 7 = log(2) for a simple Z2 spin liquid [29, 30]. 

In 2D critical systems, the behaviour of these subleading corrections to the area law 
become much more rich. One may naively suspect that, due to the diverging correlation 
length, a violation of the area law might be possible. However, it can be understood through 
a course-graining picture that the area law is still obeyed. First, assume that due to scale- 
invariance each length scale (in a renormalization group (RG) sense) contributes order 0{1) 
bit of entanglement entropy across the boundary. One may take scale-invariance to mean 
that, when rescaling the system by some factor b, the number of modes at this new length 
scale is proportional to the new boundary length i/b. Then, using this assumption and 
summing over all length scales, an entropy proportional to the boundary length is obtained. 
t 

At criticality, additional universal subleading terms to this area law are also possible, 
however they may have a complicated dependence on the geometry of the bipartition. 
Although typically believed, it is not generally known if particular geometric features, 
for example the number of vertices or the Euler characteristic [14, 31, 17, 18], give rise 
to certain universal numbers that can be compared reliably between field theories and 
quantum lattice models. This would be important in making progress towards developing 
an analog of c-theorems [32] in 2+1 space-time dimensions, which aim to identify a universal 
function with monotonic behaviour under RG fiows [33]. The fact that most field-theoretic 
calculations are limited to non-interacting systems hampers progress in this regard. In order 
to study interacting systems, one must turn to numerical techniques on finite-size lattice. 
Understandably, the geometries amenable to study on finite-sizes lattice are sometimes 
different than those that can be studied with continuum field theories, as we now discuss. 

2.1. Bipartitions with smooth boundaries 

In the continuum thermodynamic limit, dividing a systems into two partitions A and B 
is easily done with a smooth curved boundary. This geometry has become particularly 
important in entanglement monotonicity studies, where for example, in Lorentz-invariant 
theories, Casini and Huerta have shown that the entanglement of a smooth circle of 

X We thank M. Hastings for pointing out this argument, which might possibly be traced back to J. Preskih. 
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circumference ^ scales as Si ^ Al — 7', where the universal constant decreases along an RG 
flow [34]. It is thus a compelling candidate for the monotonic c-function. This may also be 
related to the entanglement of a three-sphere in odd space-time dimensions, which contains 
a subleading constant term which changes monotonically along RG flows in holographic 
theories [35, 36]. Such developments may provide a route to a (2-f-l)-dimensional analog 
to the c-theorem, especially if the fixed-point value of the monotonic quantity could be 
determined. 

In essentially any interacting theory (and some non- interacting theories) in 2+1, this 
task would necessarily fall to numerical simulations. Unfortunately, such curved geometries 
are inaccessible on lattices, where obtaining boundaries with sharp vertices is unavoidable 
when attempting to draw smooth curves. In the case where the vertices disappear in 
the continuum, it is unknown how such "pixelization" might affect the approach to the 
thermodynamic limit. In the case where vertices or corners remain in the geometry in the 
continuum, they will contribute an additional universal factor, as we discuss in Section 2.3. 

Additionally, one may also examine the entanglement entropy across a smooth boundary 
perturbed away from the critical point, where the correlation length becomes finite. In this 
case, even if the boundary has no curvature, one expects [16] , 

Sa = Ai + Ta ' 

for spatial dimensions D. Since the definition of fiat boundaries is possible in lattice systems, 
recent numerical calculations on interacting lattice models have been able to make important 
comparisons of calculated perturbatively in interacting field theories. In these numerical 
works, series expansion techniques are able to capture entanglement contributions across fiat 
boundaries on infinite lattices by systematically including larger cluster sizes [19, 13]. 

An important distinction between these numerical series expansion techniques and 
quantum Monte Carlo (QMC) is the fact that the latter is typically restricted to periodic 
finite-size systems (i.e. toroidal lattices of size L x L), meaning that these methods approach 
the thermodynamic limit in a different way. A smooth spatial boundary in a two-dimensional 
toroidal lattice is possible only if one bifurcates the torus into two separate cylinders. This 
geometry can also be studied in certain field theories, which is important for addressing the 
potential universality of cutoff-independent subleading scaling terms in the Renyi entropies. 
In this paper will we study this geometry in the context of the TFIM 2+1 QCP in great 
detail. 

2.2. Bifurcated torus: two-cylinder entropies 

As discussed above, in two spatial dimensions QMC simulations typically take place on 
toroidal lattices of size L x L. For the measurement of entanglement between two subregions 
A and B, a simple geometry is illustrated in Fig 1, where upon bifurcation two cylinders of 
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Figure 1. A toroidal simulation lattice, divided into two entangled cylinders A and B, of 
size X X L, and {L — x) x L. 

"length" X and L — x are produced. Past numerical studies of strongly-interacting gapless 
systems have demonstrated that this geometry is a sensitive probe of the entanglement 
structure of the wavefunction [37, 38, 22]. Two possible features are particularly prominent 
as one varies the length x: a striking even-odd effect which arises due to dimer-like physics in 
the wavefunction; and, a smooth a;-dependent curvature that may contain universal scaling 
behaviour. This smooth x-dependent curvature is not present in gapped states, making it a 
sensitive indicator for gapless behaviour in general wavefunctions [38] (see Section 2.2.2). 

2.2.1. Even-odd effect: As first observed in Ref. [38], a striking "even-odd" branching effect 
is observed in Sn{x/L) in certain RVB wavefunctions. In Ref. [22], this effect was understood 
in a free scalar field theory in the 2+1 dimensional quantum Lifshitz model, which is the field- 
theory of a Rohksar-Kivelson (RK) Hamiltonian. It arises due to the underlying dimerization 
of the wavefunction - and not from any underlying symmetry breaking, as one might naively 
expect for example in a valence-bond solid phase. 

Away from this exactly-soluble free field theory, the even-odd effect serves as a sensitive 
probe of the degree of dimerization in the low-energy effective description of the wavefunction. 
For example, the Neel ground state of the 2D spin-1/2 Heisenberg model can be described 
in an RVB singlet-basis where long singlets, decaying as occur [39]. Correspondingly, 
no even-odd branching effect is observed in the Renyi entropy 5*2 [37]. Similarly, for QCPs 
which are described by theories "sufficiently" far from RK-like Hamiltonians, it's reasonable 
to expect that no even-odd branching occurs. As we will see in this paper, in the 2+1 
dimensional QCP in the transverse-field Ising model (TFIM), no such even-odd branching 
occurs. 

2.2.2. Length dependence from conformal field theories: Previous numerical studies of 
interacting models on two-cylinder entropies show a clear geometry-dependent function, that 
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occurs even in the presence of even-odd branching (described above), and which depends on 
the cyhnder length x. Previous resuhs on the Neel ground state of the Heisenberg model, 
and the square-lattice RVB wavefunction showed reasonably good numerical fits [38] to the 
shape-dependence motivated from 1+1 CFT in Eq. (1), specifically, 

Sa = A^ + h log(£) + c log(sin(7r?/)) + d, (6) 

where however 6 7^ c in general. The term proportional to log(sin(7ry)) is heuristically 
included in order to account for a strong shape-dependence observed in these gapless 
wavefunctions in 2+1, in analogy to the 1+1 exact result. Importantly, a non-zero 
logarithmic term h of order unity was first observed in QMC data on the Heisenberg 
model in Ref. [37]. This phenomena was subsequently explain by Metlitski and Grover 
[40] as arising from the presence of Goldstone modes and the restoration of symmetry in a 
finite-volume simulation cell. The coefficient of this additive logarithm should be universal, 
h = Ng{D — l)/2, where Nq is the number of Goldstone modes. Thus, it is only expected to 
be present in the case of spontaneously broken continuous symmetry. It has previously been 
demonstrated not to exist in an exactly-solvable finite-size model of free spinless fermions 
on a square lattice, with n flux through each plaquette [38]. As we will see in Section 4.1, 
QMC data for the TFIM quantum critical point is also consistent with 6 = 0. 

Recently, motivated by the study of dimer RVB wavefunctions that in the continuum 
limit have conformal invariance in 2D space, CFT techniques have been used to derive an 
alternate scaling form of the shape-dependent piece for Fig. 1 [22, 41]: 

s^= Ai + cUy) + ---, (7) 
Ja{y) = z log 



1-a " L^3(2r)03(r/2) v{'2yrM2{l - y)r) J ' 

where y = x/L is the fractional width of the strip (referred to as iy/L + y in some of the 
previous literature), 1] is the Dedekind eta function and 61, is the Jacob i-Theta function. The 
above form for Ja{y) apphes for the Renyi entropies with a > 2, but the definition does not 
extend to the von Neumann entropy. Through the rest of the paper we use J{y) = J2{y) for 
simplicity. Due to our geometry (always using Lx L systems), r = iL^/Ly = i never changes 
in the above equation. Although this universal function was derived for the quantum Lifshitz 
field theory, it is interesting to test its universality away from the critical points describing RK 
and RVB-like wavefunctions. In the present paper, we explore its potential for universality 
by comparing the functional dependence of the Renyi entanglement at the TFIM QCP to 
fits of Eq. (6) and Eq. (8). 

2.3. Polygons on a torus: entanglement due to vertices 

It is well-known that, in critical systems in 2+1 dimensions, another cut-off dependent 
contribution to the entanglement entropy distinct from the boundary ("area" -law) is induced 



8 



Figure 2. A toroidal simulation lattice with a rectangular entangled region A with four 
90° vertices. 

by the presence of vertices or corners. This term can be seen to have a logarithmic dependence 
on £ through heuristic mode-counting arguments in a renormahzation group framework. 
To begin, one assumes that each length scale (in an RG sense) contributes C(l) bit of 
entanglement entropy for each vertex or corner. Then, summing the contributions from each 
length scale, one discovers that the corners contribute a constant times the total number of 
length scales, log(£). § 

Refinements and generalization of this argument to other geometrical boundaries exist 
in the literature. Most promising for comparison between analytical field theory and lattice 
numerics is the simple 90-degree vertex. In Ref. [15], the full angle-dependence of the 
logarithmic vertex contribution is explored, and found to obey, 

5« = A£ + n,a„(^)log(£) + ---, (9) 

where ric is the number of vertices, and aa{9) is a universal quantity independent of lattice 
cutoff. For the square lattices studied in this paper, it is natural to use a vertex of = 7r/2. 
The quantity aQ,(7r/2) will be universal between the continuum field theory and the lattice 
model, provided that the angle in the lattice entangled region approaches a 90° corner as 
one increases the number of sites to infinity. 

Past numerical studies have calculated aa{0) for a variety of Renyi indices a at the QCP 
of the TFIM. Of relevance to the present study, the value of a2(7r/2) has been calculated 
several times in the past literature. A value for the non-interacting fixed point in a scalar 
field theory was provided by Casini and Huerta, a2(vr/2) = —0.0064 [42]. Series expansion 
studies of the interacting 2+1 dimensional QCP of the TFIM give —0.0055(5) [19] while 
Numerical Linked-Cluster Expansion (NLCE) gives —0.0053; both results are consistent 
with each other, and lower than that result of the free field theory. Finally, previous finite- 
temperature QMC calculations on a single 36 x 36-size lattice report —0.0075(25) [20]. In 
the current paper, we aim to improve on these QMC results by employing a new projector 
method that converges the Renyi entropy 5*2 directly at T = 0, as now described. 

§ See footnote 1. 
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3. Projector Quantum Monte Carlo 

Due to a diverging correlation length, the numerical study of quantum phase transitions in 
strongly-interacting models requires extraordinary care [43]. Every effort must be taken 
to converge data on lattices of as large a size as possible, such that reliable finite-size 
extrapolations may take place. Quantum Monte Carlo (QMC) on sign-problem free models 
provide an unbiased numerical procedure to systematically approach the thermodynamic 
limit. For lattice spin models, such as the transverse-field Ising model (TFIM) studied here, 
the standard procedure used to access quantum critical behaviour involves using highly- 
efficient /inzte-temperature algorithms, such as continuous world-line [44] or Stochastic Series 
Expansion (SSE) [45, 46], operating at sufficiently low temperatures. Then, the temperature 
T is either decreased until it is converged to its T — )• behaviour for each lattice size (and 
parameter set) of interest; or, the inverse temperature (3 = 1/T is set proportional to the 
linear lattice size, (3 (x zL (where z is the dynamical scaling exponent), for each lattice size 
studied. 

Recently, a new flavour of QMC algorithm has emerged that combines the simplicity 
and efficiency of SSE QMC, with the ability to study ground state properties of model 
Hamiltonian directly at T = [39]. Such "projector" methods share many of the features of 
their T > counterparts, such as: the direct numerical coding of a D-dimensional quantum 
system to a D + 1 dimensional classical configuration, which is represented and sampled on 
a computer; and the existence of the prohibitive "sign-problem" for fermonic and frustrated 
systems. Unlike T > simulations, these methods do not operate in a D -|- 1 simulation 
cell that is periodic in the temporal direction; rather, they operate with a Hamiltonian 
on a trial state, repeatedly, projecting out the ground state, as explained below. These 
projector methods have been widely adopted for the study of T = properties of SU(2) 
(and SU(A^) [47]) invariant Hamiltonians with singlet ground states, providing a simple and 
efficient method to converge ground state properties [39, 48, 49]. In the next section, we 
describe a procedure, first developed by Sandvik [50, 51] (and recently generalized [52]), 
for the efficient simulation of the U(l) symmetric TFIM Hamiltonian, using an adapted 
projector QMC method which employs non-local cluster updates. In Section 3.2, we describe 
how this algorithm may be adapted to calculate the Renyi entanglement entropies using a 
straight-forward adaptation of the "replica" trick [53]. It is this method which allows us to 
accurately study the finite-size scaling of 5*2 at the QCP of the TFIM, and to access the 
universal sub leading quantities of interest. 

3.1. Algorithm for the transverse- field Ising model 

In 2005, Sandvik introduced a ground state projector QMC method for SU(2) quantum 
spins, using the so-called valence-bond basis [39]. At T = 0, QMC methods are tasked with 
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calculating the operator expectation value, 

{o) = ^{^'\om. (10) 

Here, one aims to use some procedure to sample , the ground state wavefunction of a 
Hamiltonian, where the normalization is Z = (ip^lilj^). The transition probabilities of the 
Metropolis algorithm are based on this overlap; the non-orthogonality of the valence-bond 
basis ensures that this is trivially non-zero. However, as we will see below, this trial state is 
not strictly required to be an overcomplete basis. In the case of the TFIM, an orthogonal 
basis can be used, provided that wavefunctions are sampled such that this overlap is 
non-zero, using a non-local cluster algorithm. 

In a projector QMC representation, the ground state wavefunction is estimated by a 
procedure where a large power of the Hamiltonian is applied to a trial state, call it To 
see the projection of the ground state wavefunction, one can write the trial state in terms 
of energy eigenstates {ip^), n = 0, 1, 2 . . ., |a) = ^„ Cn]'?/'"'), so that a large power of the 
Hamiltonian will project out the ground state. 



^ol 



-H)"'\a) = co\E, 

Co|-Eo|"'|V'°) as m oo. 



'IV. 



Here, we have assumed that the magnitude of the lowest eigenvalue \Eo\ is largest of all the 
eigenvalues. To achieve this, one may be forced to add a sufficiently large negative constant 
to the overall Hamiltonian (that we have not explicitly included). Then, from this expression, 
one can write the normalization of the ground state wavefunction, Z = {ip^lip^) with two 
projected states (bra and ket) as, 

Z = (a|(-/7)™(-/f)"|a) = (a|(-if)2-|a), (12) 

for large m. In a procedure that will be familiar to any SSE aficionado [54], the Hamiltonian 
is written as a (negative) sum of elementary lattice interactions 

t a 

the indices t and a referring to the operator "types" and lattice "units" over which the terms 
will be sampled. In order to represent the normalization as a sum of positive-definite weights 
we can insert a complete resolution of the identity between each Ht^^a^ 



2m 

^=EEn<"'i^*.">^>- (14) 

{a} Sm j=l 

Where this equation has been cast in a form similar to that for finite-T SSE. Note that, the 
sum over the set {a} and the operator list Sm = [ti, di], [^2, 0,2], - ■ ■ [t2m, 02m] must be done 
with importance sampling. As we will see below, an update procedure can be constructed 
that efficiently samples both the list of operators Ht^a, and (separately) the left and right 
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basis states. Thus, for sufficiently large m, any trial state \a) can be chosen, which is 
equivalent to using the equal superposition of all spin states [51]. Other choices, such as 
a variationally optimized state, may also be used [51]. 

Turning to the Hamiltonian Eq. (2), a convenient definition of operator types is, 

H^,,a = h{a+ + a;), (15) 
Ho,a =h, (16) 
Hi,a =J{atcr] + l). (17) 

Note that the index a has two different meanings: a site or a bond, depending on the operator 
type. It is evident that some simple constants have been added to the Hamiltonian (Eq. (2)): 
the diagonal operator ifo,a, and also the +1 in Eq. (17). The ffist results in matrix elements 
with equal weight for both one-site operators. Denoting each matrix element in the standard 
basis: the erf = +1 eigenstate is | • )j and cr| = — 1 is | o )j, the non-zero matrix elements 
are, 

( . o ) = ( o . ) = /i, (18) 

{.\Ho,a\*) = {o\Ho,a\o) = h, (19) 

( • • \H,^,\ • • ) = ( o o \H,^,\ o o ) = 2J. (20) 

The -D + 1 dimensional projected simulation cell is built such that 2m operators of the type 
15 to 17 are sampled between the "end points" (i.e. the trial states). Then, sampling occurs 
via two separate procedures, as follows. 

First, the diagonal update where one traverses the list of all 2m operators in sequence, 
e.g. from \ai) to ja^). If an off-diagonal operator H_i a is encountered, the spin associated 
with that site is flipped but no operator change is made. If a diagonal operator is encountered, 
the Metropolis procedure is: 

(i) The present diagonal operator, i^o.a or Hi a, is removed. 

(ii) A new operator type is chosen, t = or t = 1, corresponding to the insertion of either a 
diagonal h or a diagonal J operator, i.e. Eq.(19) or (20). The transition probability to 
add Ho,a is, 

hN 

^^^''''^ = hN + {2J)N,' ^^^^ 
where and A^^ are the number of sites and bonds in the lattice, respectively. Note, 

(iii) If ifo,a is chosen, a site a is chosen at random, and the operator is placed there. 

(iv) If Hi a is chosen, a random bond a is chosen. The configurations of the two spins on this 
bond must be parallel for the matrix element to be nonzero. If they are not, then the 
insertion is rejected. Steps (ii) to (iv) are repeated until a successful insertion is made. 
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Figure 3. Two representations of a 1+1 dimensional simulation cell with m = 6, before 
(left) and after (right) a cluster move that updates site operators and basis states. The 
basis states tr^ = 1 (cr^ = ~1) are filled (open) circles, and solid (dashed) lines connecting 
operators. Off-diagonal site operators, Eq. (18), are represented by filled squares, whereas 
diagonal site operators, Eq. (19) are open squares. Vertical bars represent bond operators, 
Eq. (20). The mid-point of the 1+1 simulation cell, discussed in the text, is indicated by a 
vertical dashed line. 



One can see that this diagonal update is necessary in order to change the topology of the 
operator sequence in the simulation cell. However, in order to get fully ergodic sampling of 
the TFIM Hamiltonian operators, one must employ non-local updates in addition to these 
simple diagonal updates. 

For T = projector methods, non-local updates have been discussed previously in the 
literature in the context of the valence-bond basis [49]. In the present case, the TFIM 
Hamiltonian does not conserve cr^, and a different type of branching non-local update, called 
a "cluster" update, must be used. We use cluster updates adapted from the finite-T SSE 
procedure described in Ref. [50]. A crucial observation is the judicious choice of -ff_i,a and 
Ho^a to both have the weight h, which allows for unrestricted sampling between the two 
operator types. One can easily see that a functional definition of non-local clusters will be 
an Ising spin forming a space-time group, bounded by either single-site operators, or by spin 
states of the end point trial states \ai) and la^) (see Fig. 3). If all spins within a cluster are 
flipped, the total weight of the configuration remains unchanged. One is then free to build all 
clusters deterministically, and flip each with a Swendsen-Wang algorithm, i.e. a probability 
of 1/2. In this way, one sees how a fully ergodic sampling of both the operator types and 
basis state is sampled in the projector QMC. 

In Fig. 3 special care has been taken to note the mid-point of the simulation projection, 
since operator expectation values are measured there: 

_ (^,^1^.°) _ {ai\{-HrO{-Hr\ar) 

It is particularly important to reiterate, especially in anticipation of the next section, that 
the wavefunction overlap (V'^'lV'r) must be non-zero (indeed, the Metropolis sampling is set 
up to force this). This is ensured by requiring that spin states are the same on the left and 
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Figure 4. An illustration of the "SWAPa" operator acting on a replicated simulation cell, 
for the calculation of 82- Basis states and TFIM operators follow the same legend as Fig. 3. 
Here, each simulation has three physical spins, and m — S operators in each non-interacting 
replica, represented one on top of the other. At left, clusters within each replica that cross 
the mid-point of the 1-1-1 simulation cell are highlighted in different colours. At right, a 
SWAPyi operator has been applied to permute the basis states corresponding to region A, 
between each replica. The reconfiguration of the clusters crossing the mid-point can be seen 
from their colours. 



the right of the "mid-point" of the D + 1 dimensional celL Since spins are only constrained 
to be the same within a given cluster, a proper normalization for expectation values could 
be defined as, 

Z = = 2^0, (23) 

where A'^o is the number of independent clusters that cross the boundary. This follows from 
the fact that each connected cluster in Fig. 3 has two possible spin orientations, = 1 and 
= —1, independently. 

Various relevant measurements, such as the ground state energy, can be constructed for 
this simulation [51]. However, for the purposes of this paper, we are interested specifically 
in the Renyi entanglement entropies, which only require knowledge of the cluster structure 
of the simulation at the middle of the propagated simulation cell; albeit in a non-trivial 
replicated (or multi-sheet) geometry. Thus, in the next section, we describe the specific 
implementation of the Renyi entropy estimator for the projector QMC for the TFIM. 



3.2. Measuring Renyi entropies through the SWAP a operator 



In a development crucial to the study of entanglement at interacting quantum critical points, 
QMC methods have recently been introduced that are capable of measuring the degree of 
entanglement in a ground state wavefunction, specifically using the Renyi entropies [1], 

1 



Sa 



1 — a 



log 



Tr(p; 



(24) 



14 



for integer a>2. Here, is the reduced density matrix, = Tr^jp}, and A is a subregion 
of a lattice system (with B being its complement) - see Figs. 1 and 2. The von Neumann 
entanglement entropy corresponds to the limit a — )■ 1. 

Unlike a linked-cluster expansion or other methods based on Lanczos diagonalization 
[13], the direct measure of Renyi entropies can not be done in QMC using conventional 
estimators. However, recent work has demonstrated that it is possible to measure for 
integer a > 2 using a replica trick [6, 14, 55, 56, 57]. For T = 0, the calculation of Sa is done 
via the following procedure: 

(i) The system is copied into a independent "replicas" . 

(ii) Each replica is independently projected to sample the ground state. The total 
wavefunction for the system of all replicas is denoted by 

(iii) Following Eq. (22), the operator O is replaced by the "SWAP^" operator for a = 2 
[53] or permutation operator for a > 3 [37]. This operator literally swaps (or permutes) 
basis states in the region A between the a copies. See Fig. 4 

(iv) The expectation value 

(SWAP,) = Wl^^'^f^l'^') = 2"--«, (26) 

where Na is the number of independent clusters crossing the middle of the re-connected 
( "swapped" ) partition function, and Nq is the number of clusters before the swap took 
place. Note: here Nq is the number of clusters in all a replicas or copies. 

(v) Finally, Eq. (24) is applied: e.g. ^2 = - log [(SWAP^) 

Note that. Step (iv) comes from the fact that Eq. (22) is used, with O = SWAP^ as 
the operator being measured. This results in the simple procedure of measuring the overlap 
in the numerator and denominator - both of which are simply the number of spin states per 
cluster (2) raised to the power of the number of independent clusters crossing the middle of 
the -D + 1 dimensional simulation cell. 

As first pointed out in Ref. [53], as the lattice size (and particularly the size of region A) 
grows, sampling statistics become exponentially poor when using the naive SWAP^ operator 
as described above. Hence, a slight adaptation called the ratio trick must be used in order 
to improve statistics. The ratio trick involves calculating the Renyi entropy in several steps, 
each step being a separate simulation that involves sampling a ratio: 
(V.°|SWAP^|V;°) ^ 

{^flSWAFxm ^ ' 

where X denotes a subregion that is spatially smaller than A. In other words, the weight 
of a simulation becomes (un-physically) related to a partially-swapped simulation cell. In 
a practical simulation, many spatial sub-regions Xi, each employed as the weight of a 
separate simulation, are used to build up towards the physical bi-partition A of interest. 
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The Renyi entropy is then built up by multiplying the contribution of Eq. (26) from each 
spatial subregion. 

For the TFIM simulation discussed in this paper, the procedure for efficiently calculating 
the Renyi entropy closely follows that used in another T = projector QMC - the valence- 
bond basis QMC for the spin- 1/2 Heisenberg model [39, 53], which has a detailed description 
(including the ratio trick) in Ref. [37]. Remarkably, the only algorithmic difference between 
the two QMC algorithms is the structure of the space-time clusters employed in the projector 
method. In Eqs. 25 and 26, the A^-numbers {N^, Nq, Nx) count the number of branching 
clusters, spanning both the spatial and propagation directions, that cross the middle of the 
simulation cell. In the spin-1/2 Heisenberg model [37], this number counts the non-branching 
loop structures that cross this middle point, due to the different nature of Hamiltonian 
operators in that model. As in the case of the TFIM clusters, each Heisenberg loop in 
the valence-bond representation has two spin states associated with it (for SU(2); this is 
modified to for SU(A^)). It is remarkable that Hamiltonians with different symmetries, 
and completely different basis-state representations in the projector QMC, end up with 
equivalent measurement procedures for the Renyi entanglement entropies. 

3. 3. Convergence of S2 at a quantum critical point 

In this paper, we examine the second Renyi entropy, setting a = 2 in Eq. (24). Like any 
other observable measured with the T = projector method, the value of 5*2 will have its 
own unique convergence properties, as a function of the projector length m, for each value 
of h/J studied. When focussing on the critical point, h/J = 3.044, where the correlation 
length diverges, one must be particularly careful to ensure that the simulation is converged 
in m for each lattice of size N = L x L. 

In Fig 5, we examine the value of 5*2, using the two-cylinder geometry of Fig. 1, at the 
point when the two entangled cylinders are of equal size, x = L/2. Of the geometries studied 
in this paper, this point is expected to be the most difficult to get good statistics on, and we 
use the ratio trick, Eq. (26) to converge it, building up each subregion Xj with at most 1 x L 
sites (less for larger sizes). See also the discussion in the next paragraph for a comparison 
to data taken without the ratio trick. For each of the four system sizes that we have studied 
in detail, we see that the value of 5*2 at the centre-point x = L/2 converges for sufficiently 
large m - requiring a slightly larger value of m/N as A^ increases. For the largest system 
size that we collect detailed convergence data on, L = 22, the value of S2 saturates between 
50000 < m < 60000 operators, i.e. m/N slightly larger than 100. We continue to collect 
data (see Fig 5) for larger m, but find that the value of S2 is essentially converged within 
error bars. For the data collection in the rest of the paper, we settled on a fixed m/N = 160, 
which is more than sufficient to converge L = 22 and larger. A comprehensive study of the 
m-dependence of the Renyi entropy for different values of a and A^ would be an interesting 
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Figure 5. The convergence of the second Renyi Entropy, S2{x — L/2), at the centre point 
of the two-cylinder geometry compared to the value at the maximum m tested, mniaxi as a 
function of the number of sites over operators applied to the trial state, N/m. Simulations 
are performed at the dashed line, m/N = 160. Inset: The second Renyi entropy versus the 
cylinder length, for = 24 x 24. Note the significant increase in the statistical quality of 
the sampling when the ratio trick, Eq. (26), is employed. 



topic of future study. 

In the inset of Fig. 5, we see raw data for the x-dependence of 5*2 using the geometry in 
Fig. 1, on an = 24 X 24 lattice with m/N = 160. There, we have compared data obtained 
from a single simulation using a bare measure of the SWAP^ operator, Eq. (25), with that 
obtained from a procedure where the ratio trick, Eq. (26), is used to build up each entangled 
region A. Clearly, naive measurement of the bare SWAP^ operator is insufficient to get 
controlled statistical sampling for large x values. We find that, in the case of the TFIM at 
h/J = 3.044, the ratio trick is absolutely necessary for simulation of size L > 16. 

4. Simulation results on finite-size lattices 

In this section, we report results for the second Renyi entropy, 5*2, for the Hamiltonian 
Eq. (2) precisely at the QCP {h/J = 3.044) using the QMC simulation method discussed 
in the last section. We discuss two entangling geometries: a bipartition between A and B 
that smoothly cuts each torus into two cyhnders (Fig. 1), and a square bipartition with four 
90-degree corners (Fig. 2). 
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Figure 6. (Color online) The residual of fitting all of the x — L/2 and a; = L/4 strips to 
the form 5*2 (i) = aL -t- d, keeping in mind £ — 2L for both geometries. From this we see 
that the data fit is consistent with Eq. (6), with b — (i.e. no log(^) dependence). 



4.1. Bifurcated torus: two- cylinder entropies 

The first geometry we consider is that of an L x L torus, where the two entangled regions 
result from smoothly cutting the torus into two cylinders, as shown in Fig. 1, where the 
length of each cylinder is x x L and {L — x) x L. 

Before examining the full x/L dependence of the two-cylinder geometry, we discuss the 
possibility of a non-zero b in Eq. (6) by examining regions A of a fixed x/L embedded in 
different finite-size lattices L. In Eq. (6), one may eliminate the x-dependence of the term 
proportional to c by fixing x/L to be a constant; e.g. if x is L/2 or L/A, this contribution will 
be absorbed into the additive constant d. Fig. 6 shows the residuals, that is 5*2 (L) — aL — d, 
for two different choices of x/L as a function of system size for the half and quarter cylinder 
partitions, with b explicitly set to zero. From this, we conclude that the fit is acceptable 
within statistical errors. If, instead, we allowed a 6 7^ as a fit parameter, a very small 
negative coefficient is found, b ^ —0.01, two orders of magnitude smaller than that found 
in systems with continuous symmetry breaking [37]. As such, we argue that this value is 
inconsistent with both physical expectations (see Sec. 2.2) and the results below (see Fig. 8). 
We conclude that our data is consistent with the absence of a subleading log{£) dependence 
in the case of smooth boundaries. 

In the following, we assume that no additive logarithms are present in the entanglement 
scaling at the TFIM quantum critical point; 6 = in Eq. (6). Fig. 7 then shows the 
entanglement entropy as a function of the cylinder size, x, for our three largest system 
sizes, L = 28, 32 and 36. Recall that in some other gapless wavefunctions, e.g. the RVB 
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wavefunction studied in Refs. [38] and [22], a very prominent branching effect is apparent 
whicli produces separate entanglement curves for even-x and odd-x (see Sec. 2.2.1). As 
discussed by Steplian et al. [22], tliis even-odd effect is a measure of liow RVB-like a 
wavefunction is. In the present case, it is clear that if any even-odd effect occurs, it is 
beyond our ability to detect in these QMC simulations. 

Next, we are interested in examining the functional dependence of the shape of the 
curves in Fig. 7. To do this, we look at a variety of system sizes, L, and a variety of 
entangled-region lengths, x, and examine the fit of the entropy to two equations. The first is 
the form log(sin(7rx/L)), Eq. (1), which is relevant for 1+1 dimensional systems, and which 
has been used in an ad hoc way in the past analyze some 2+1 dimensional entanglement 
entropy data [21] including systems with a 2+1 dimensional fermi surface [38] where Eq. (1) 
gives a reasonable (but not perfect) approximation to the fit. The second is the RVB shape 
function derived by Stephan et al. [22], Eq. (8). For the purpose of a numerical comparison, 
we fit to functions of the form, 

S'log = ai + CL log(sin(7ry)) + d, (27) 
S^yB = a^ + CLJ{y) + d, (28) 

with y = x/ L and where the constants a, cl, and d may be different for the above equations. 
Note that we allow the freedom that cl can vary with system size. However, a and d are fit 
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Figure 8. (left) The residual of fitting the entanglement entropy to (orange) Eq. (27) 
and (teal) Eq. (28) for all system sizes. Note the systematic deviation of the redisual in the 
case of the log fit, that is, the strong tendency to over estimate (negative residual) the fit 
at small x/L, and under estimate (positive residual) the fit at x/L near 0.5. (right) The 
minimizing parameter cl for the two fits, with a non-trivial size dependence in the case of 
the log fit. 

by considering all systems sizes together. In this way, we can use the variation of cl with L 
to examine the quality of fit for each of the two functions. 

Fig. 8 shows the residual, S{L, x) — S\og and S'(L, x) — Srvb, of the fit of the entanglement 
entropy to Eqs. (27) and (28). In the right panel a comparison of the minimizing cl for both 
functional fits is shown. As is evident from the left-hand plot, the residuals are comparable 
for the two functional forms, with Eq. (28) giving a slightly better fit, using this metric. 
However, an important point to notice on the right-hand plot is that, when we attempt 
to fit the data for all system sizes L to the form Eq. (27), there is a clear system-size 
dependence in c^. In Ref. [21], using geometries amenable to DMRG simulation, it was 
argued that this increasing trend of cl may level off at a moderate system size value; this 
appears to support the notion of a central charge c = 1 in the limit of large system size. As 
is evident in Fig. 8, no such conclusion can be drawn from the toroidal QMC geometries for 
the second Renyi entropy - although, the QMC is not able to probe the behaviour for the 
von Neumann entropy, making a direct comparison difficult. In QMC there is no evidence 
that Cl is converging to a constant value for Eq. (27), leading one to instead conclude that 
this cannot be a system-size independent (and hence universal) quantity in 2-1-1 dimensions 
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Figure 9. The scaling of an £/2 x £/2 square (of boundary length £) in a 36 x 36 system. 
The data is subtracted in such a way that only the log term should remain, with the line 
showing the best fit to this assumption. Note that the best fit of the entire data set gives 
a positive a2(7r/2) and there is systematic curvature away from a log, that is to say the 
deviation of the data seems correlated in x rather than random and gaussian. i here is 
defined as the (average of inner and outer) boundary length dividing the two regions. 

for the Renyi entropy. 

In contrast, for the case of Eq. (28) there is a much less-pronounced size dependence for 
the coefficient of J{y). As evident in Fig. 8, a much better case can be made that cl levels 
off to a system-size (ie. cutoff) independent value in the limit of large lattice sizes. This fuels 
speculation that J{y) may be a universal scaling function relevant for all fixed points in 2+1 
dimensions, as discussed more in Section 5. 

4-2. Polygons on a torus: entanglement due to vertices 

The second geometry we examine is that of a rectangle embedded in a torus, as shown in 
Fig. 2. As discussed in Section 2.3, this geometry is particularly promising for accessing 
universal subleading corrections to the area law that may be computed in both continuum 
field theories and lattice models as in this work. Since these universal corrections arise due 
to vertices (or corners) in the entangled region, any geometry (or aspect ratio) of rectangle 
should have four separate vertex contributions, 4a2(7r/2), from Eq. (9), assuming that each 
rectangle is large enough that interactions between corners may be neglected. In our QMC 
simulation geometries, the use of differently-shaped rectangles gives us more than one avenue 
to extract the coefficient of the subleading log, 02(71/2), thus providing an independent test 
for universality. 
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Figure 10. The comparison of scaling an L/2 x L/2 square for two definitions of the 
boundary length. "Square" assumes the boundary is the average of the number of sites on 
the inside and the outside of the boundary, or simply the length of each edge times four. 
"Small Square" only counts the number of unique sites bordering in the inside edge of the 
boundary, which is smaller than the previous definition by four sites. By excluding smaller 
system sizes the extracted coefficient using these two methods converges, but the error in 
the extraction also increases. 



There are two scalings that we test: Fig. 9 shows a square scaled in a fixed size 
simulation, while Fig. 10 shows an L/2 x L/2 rectangle scaled in an L x L system. Scaling 
ofanL/2xL/4 rectangle was also examined, and results are consistent with the L/2 x L/2 
rectangle within error. We generate the data using the different geometries, then fit it to 
Eq. (9) with an additional allowance for a constant term in the fit: 



Repeating a previous analysis [20], we first perform a preliminary fit using a fixed lattice 
size, L = 36, and examine the Renyi entropy of all possible squares up to size L/2 x L/2. 
This data is fit to the form Eq. (9), and the scaling of the subleading logarithmic term as a 
function of boundary length is extracted. The result of this analysis is plotted in Fig. 9; we 
find that the constant is actually of the opposite sign from that seen in previous work [20] 
(refiecting a significant decrease in the error bars on S2 with the present data set). If we 
exclude smaller rectangles from the fit, the sign of 02(77/2) remains opposite to previous 
work. In addition, if we look at the quality of the fit in Fig. 9, there is a large systematic 
deviation from a log{£) form, suggesting this geometry may be infiuenced by other factors 
not accounted for by this analysis, in particular, terms which depend on the aspect ratio of 
region A. 

In Fig. 10, we attempt another type of analysis aimed at eliminating the systematic 
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shape-dependence observed in Fig. 9. There, the size of region A is fixed to be a square or 
rectangle with a perimeter proportional to the system's size L x L. Then, many different 
LxL toroidal lattice sizes are individually simulated (L = 10, 12, 14, 16, 20, 24, 28, 32, 36, 40). 
Fig. 10, which shows the subleading logarithm dependence isolated from the area law (and 
additive constant), demonstrates that this analysis produces a very good fit to the expected 
function log(£). Analyses of two geometries {L/2 and L/A) of region A produce a consistent 
value for the universal coefficient of this logarithm, with the value extracted for squares of 
size L/2 X L/2 being more accurate (due to the availability of more system sizes). 

There are two additional pieces of analysis required for completeness: size exclusion and 
the definition of i. For all results presented thus far on the polygonal entangling geometry, 
we have taken the definition of i to be the average of the number of sites on the inside 
and outside of the boundary, except where specifically indicated otherwise. An alternate 
definition of the boundary length uses the minimum of the number of sites on the inside 
and outside of the boundary. These definitions of the boundary length only differ when 
considering our rectangular regions, differing by four lattice spacings (assuming a square of 
at least 2 x 2 in size). With very large sets of data, size exclusion can be used to eliminate bias 
when fitting data for which subleading corrections are not known (here, those proportional 
to 1/L). Since the number of available sizes in our study is limited, we can only do a limited 
exclusion study. The result of this analysis is that the corner term tends to become larger 
as smaller systems are excluded, and with smaller systems included the smaller definition of 
i suggests a smaller (in magnitude) value of a2(vr/2). It should also be noted that with a 
less comprehensive analysis, it is possible to get values with an artificially lower error bound, 
and part of the reason for this work is to illuminate the possible pitfalls in the extraction of 
these universal subleading terms. 

All of this analysis suggests a universal corner contribution of a2{TT /2) = —0.006(2), with 
caveats mentioned in the previous paragraph, for the TFIM at criticality, i.e. the universality 
class of the Wilson-Fisher fixed point. Remarkably, this value is very close to the value 
calculated in a continuum field theory at the non-interacting Gaussian fixed point by Casini 
and Huerta, a2(vr/2) = —0.0064 [15, 42]. Previous numerical series expansion studies of 
the interacting 2-1-1 dimensional QCP of the TFIM give —0.0055(5) [19] while Numerical 
Linked-Cluster Expansion (NLCE) gives —0.0053 [13]; both results are consistent, and lower 
than the result from the free field theory. Our QMC result gives independent confirmation 
of the vahdity of the series and NLCE results; however, due to statistical error bars, it is 
impossible for the QMC to distinguish between the non-interacting fixed-point value and the 
previous estimates for the interacting theory. 
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5. Discussion 



Using a novel "projector" quantum Monte Carlo (QMC) method that operates at zero 
temperature [51], we have performed a detailed numerical study of the Renyi entanglement 
entropy, 5*2, at the quantum critical point of the transverse-field Ising model in 2+1 space- 
time dimensions. We have focussed on two different entangling geometries amenable to 
large-scale simulation on toroidal lattices of size L x L. 

First, when the entangling region is a square or rectangular polygon with four 90-degree 
corners, we confirm the expected scaling form of 

S2 = Ae + 4a2 (7r/2) log(£) + d, (30) 

where A is the non-universal coefficient to the area (or boundary) law, and, 

a2(7r/2) = -0.006(2), (31) 

is our best estimate of the universal coefficient of the subleading logarithmic term, which 
arises from each corner. It is a very non-trivial success that this value is within error bars of 
the value calculated at the same quantum critical point by numerical series and linked-cluster 
expansions [19, 13], both of which take very different approaches to the thermodynamic limit. 
This value of the universal coefficient is also very close to the value calculated in a continuum 
field theory at the non-interacting Gaussian fixed point, a2(7r/2) = —0.0064 [42]. The reason 
for this coincidence may be the fact that the Wilson-Fisher fixed point describing criticality 
in the TFIM may be reached perturbatively from the Gaussian fixed point. We hope that 
this fact motivates future analytical calculations of aa{6) in continuum field theory at the 
interacting Wilson-Fisher fixed point. 

Second, we consider the lattice divided into two cylinders of size x x L and {L — x) x L, 
separated by two smooth (vertex-less) boundaries of length L. With this geometry, our data 
is consistent with the absence of any "additive logarithm", i.e. no logarithmic divergence 
depending on L/a. This is the same conclusion found in a previous calculation of free 
spinless fermions on finite-size square lattices, with vr-flux through each plaquette [38]. In 
addition, no even-odd branching effect, like that observed for RVB-like phases [38, 22], is seen 
in this geometry at the TFIM quantum critical point. Instead, we find excellent functional 
fit to 

S2 = Ae + cJ{y) + d, (32) 

where y = x/L is the aspect ratio of a cylinder, and J{y) (Eq. 8) is a function first 
derived for the quantum Lifshitz fixed point which describes certain Rokhsar-Kivelson (RK) 
Hamiltonians. We argue that this function gives qualitatively better fits to our QMC data, 
and additionally is consistent with a system-size independent coefficient c, which is not 
the conclusion one can draw if the data is fit to the familiar form derived for 1+1 CFTs, 
log(sin(7r?/)). It is interesting to speculate that this new function, J{y), may be a universal 
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scaling function relevant for all fixed points in 2+1 dimensions, a conjecture that could be 
addressed with other numerical and field-theoretic studies. In particular, QMC simulations 
should have the ability to calculate the numerical value of the coefficient c at different 
interacting quantum critical points, as demonstrated by the current study. Finally, we stress 
the fact that the RVB-shape function J{y) is only applicable for Renyi entropies of order 
a >2, i.e. not the von Neumann entropy, thus emphasizing the practical utility of S2 in the 
study of universality at quantum critical points. 

We hope that this work, which was obtained with a modest amount of CPU time (fa 300 
CPU-years) , will motivate the calculation of similar quantities related to Renyi entanglement 
entropy at a variety of critical points in 2-1-1 and higher dimensions in the near future. In 
addition, we hope that this demonstration of the most convenient geometries for the study 
of entanglement in finite-size lattices via quantum Monte Carlo simulations will lead to field- 
theoretical studies of similar entanglement quantities at a variety of fixed points. If such a 
synergy could be accomplished between analytical theory and numerical simulation, the past 
successes in studying universality through entanglement in 1-1-1 may soon be translated to 
2-1-1 and higher-dimensional quantum critical points. 
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